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Abstract 

Electrical activity plays a pivotal role in glucose-stimulated insulin secretion from pancreatic /3-cells. Recent findings have 
shown that the electrophysiological characteristics of human /i-cells differ from their rodent counterparts. We show that the 
electrophysiological responses in human /J-cells to a range of ion channels antagonists are heterogeneous. In some cells, 
inhibition of small-conductance potassium currents has no effect on action potential firing, while it increases the firing 
frequency dramatically in other cells. Sodium channel block can sometimes reduce action potential amplitude, sometimes 
abolish electrical activity, and in some cells even change spiking electrical activity to rapid bursting. We show that, in 
contrast to L-type Ca 2+ -channels, P/Q-type Ca 2+ -currents are not necessary for action potential generation, and, 
surprisingly, a P/Q-type Ca 2+ -channel antagonist even accelerates action potential firing. By including SK-channels and 
Ca 2+ dynamics in a previous mathematical model of electrical activity in human /i-cells, we investigate the heterogeneous 
and nonintuitive electrophysiological responses to ion channel antagonists, and use our findings to obtain insight in 
previously published insulin secretion measurements. Using our model we also study paracrine signals, and simulate slow 
oscillations by adding a glycolytic oscillatory component to the electrophysiological model. The heterogenous 
electrophysiological responses in human /i-cells must be taken into account for a deeper understanding of the 
mechanisms underlying insulin secretion in health and disease, and as shown here, the interdisciplinary combination of 
experiments and modeling increases our understanding of human [j-ce\\ physiology. 
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Introduction 

Glucose-stimulated insulin secretion from human pancreatic /?- 
cells relies on the same major signaling cascade as their rodent 
counterparts, with electrical activity playing a pivotal role. 
Following metabolism of the sugar, ATP-sensitive potassium 
channels (K(ATP)-channels) close in response to the elevated 
ATP/ADP-ratio, which triggers action potential firing and Ca 2+ 
-influx through voltage-gated calcium channels. The resulting 
increase in intracellular calcium leads to insulin release by Ca 2+ 
-dependent exocytosis [1—4]. However, the electrophysiological 
properties of human and rodent /?-cells show important differenc- 
es, e.g., with respect to their palette of expressed Ca 2+ -channels 
and the role of Na + -channels, which contribute to electrical 
activity in human but not in rodent /3-cells [1,3]. 

Mathematical modeling has played important roles in studying 
the dynamics of electrical activity in rodent /?-cells [5,6], and could 
plausibly aid in understanding the electrophysiological responses 
in human /J-cells, and how they might differ from rodent cells. 
Recently, the first model of electrical activity in human /?-cells [7] 
was constructed from careful biophysical characterizations of ion 
channels in human /?-cells, mainly from Braun et al. [3], The 



model [7] included Na + -channels, three types of Ca 2+ -channels, 
an unspecified leak-current, and several K + -channels: delayed 
rectifier (Kv) K + -channels, large-conductance (BK) Ca 2+ 
-sensitive K + -channels, human ether-a-go-go (HERG) K 
-channels as well as K(ATP)-channels. Recendy evidence for 
small conductance (SK) Ca 2+ -sensitive K + -channels in human 
/?-cells was published [4,8], a current not included in the 
mathematical model [7]. 

The model [7] was shown to reproduce, depending on 
parameter values, spiking or rapid bursting electrical activity, 
which could be modified in accordance with a series of 
experiments by simulating pharmacological interventions such as 
ion channel blocking. These experiments were in general 
straightforward to interpret, also without a model. For example, 
the facts that blocking depolarizing Na + - or Ca 2+ -currents 
slowed or abolished electrical activity [3] are as one would expect. 

Here, we extend the previous model for human /?-cells [7] by 
including SK-channels and Ca 2+ dynamics, and show that the 
model now has reached a level of maturity that allows us to get 
insight in less intuitive experimental findings. We find experimen- 
tally that SK-channels in some cells play an important role in 
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Author Summary 

Insulin is a glucose-lowering hormone secreted from the 
pancreatic /J-cells in response to raised plasma glucose 
levels, and it is now well-established that defective insulin 
secretion plays a pivotal role in the development of 
diabetes. The /3-cells are electrically active, and use 
electrical activity to transduce an increase in glucose 
metabolism to calcium influx, which triggers insulin 
release. Experimental and theoretical studies on /J-cells 
from rodents have provided valuable insight in their 
electrophysiology. However, human /3-cells differ from 
their rodent counterparts in several aspects including their 
electrophysiological characteristics. We show that the 
electrophysiological responses in human /J-cells to a range 
of experimental manipulations are heterogeneous. We 
extend a previous mathematical model of electrical activity 
in human /J-cells to investigate such heterogeneous and 
nonintuitive electrophysiological responses, and use our 
findings to obtain insight in previously published insulin 
secretion measurements. By adding a glycolytic compo- 
nent to the electrophysiological model, we show that 
oscillations in glucose metabolism might underlie slow 
oscillations in electrical activity, calcium levels and insulin 
secretion observed experimentally. We conclude that the 
interdisciplinary combination of experiments and model- 
ing increases our understanding of human /i-cell physiol- 
ogy and provides new insight in /3-cell heterogeneity. 

controlling electrical activity, while they have virtually no effect in 
other cells. Using the extended version of the model, we show that 
this difference can be explained by differences in the excitability of 
the cells. Moreover we find that SK-channels can substitute for 
HERG-channels in controlling rapid bursting. We also show that 
blocking Na + -channels in some cells can transform spiking 
behavior into rapid bursting, in contrast to the usual effect of Na + 
-channel blockers, which in general reduce or abolish spiking 
behavior [3,9,10]. Using our model we suggest that this happens in 
cells with a large Na + -current and that BK-channels play a 
prominent role. In addition, we suggest that SK-channels might 
underlie the surprising result that blocking depolarizing P/ Qjtype 
Ca 2 + -channels enhances electrical activity, in contrast to the effect 
of L- or T-type Ca 2+ -channel antagonists, which reduce 
excitability and electrical activity [3]. Our model is then used to 
investigate paracrine effects of y-aminobutyric acid (GABA) and 
muscarinic signaling on electrical activity. Finally, we show 
experimentally slow oscillations in electrical activity that might 
underlie pulsatile insulin secretion from human pancreatic islets, 
and by adding an oscillatory glycolytic component [1 1] to the 
electrophysiological model, we simulate such slow bursting 
patterns. 

Results 

To investigate a series of experimental observations, we have 
extended our previous model of electric activity in human /f-cells 
[7] by including several additional components of human /J-cell 
physiology, as described in the following, and in greater details in 
the Methods section. The mathematical modeling was carefully 
based on experimental data, as was the development of the core 
electrophysiological part modeled previously [7]. The extended 
model includes small conductance Ca 2+ -activated K + -channels 
(SK-channels), which are expressed in human /J-cells [4,8]. The 
size of the SK-current was estimated from experimental measures 
[8] . We made a special effort to carefully model the submembrane 



dynamics of Ca 2 + , since SK-channels are controlled by the 
submembrane Ca 2+ concentration ([Ca 2+ ] mem ), which reacts 
rapidly to each action potential so that activation of SK-channels 
might influence the generation and shape of action potentials 
during spiking electrical activity. In order to study paracrine 
signalling, our extended model also includes currents due to y- 
aminobutyric acid (GABA) activation of GABA^ receptors, which 
are ligand-gated CU channels operating in human /J-cells [12]. 
Finally, a glycolytic oscillator [1 1] has been added to the model to 
account for slow oscillations in ATP levels in human /J-cells 
[13,14], which have been suggested to underlie slow patterns of 
electrical activity, Ca 2+ oscillations, and pulsatile insulin release. 

Summarizing, the new version of the model now includes 
components from glucose metabolism, additional electrophysio- 
logical components (SK-channels and GABA^ receptors), and 
Ca 2+ dynamics, leading to a global model of human /?-cell 
physiology, which, importantly, is based as far as possible on 
published data from human /J-cells. 

SK channels 

When stimulated by glucose, human /J-cells show electrical 
activity [1,3]. Human /J-cells express SK-channels [4,8], which 
might participate in controlling electrical activity. To study the 
role of SK-channels in human /J-cells, we included SK-channels 
and Ca 2+ dynamics in our previous model [7]. The new model 
with standard parameters produces spiking electrical activity 
(Fig. 1A), which is virtually unaffected by setting the SK- 
conductance gsK = 0 nS/pF simulating SK-channels block. This 
model prediction was confirmed by our experimental data, and 
was also observed in at least one cell byjacobson et al. [8]. Fig. IB 
shows an example of spiking electrical activity in a human /J-cell 
stimulated by 6 mM glucose, where addition of the SK1-3 channel 
blocker UCL 1684 (0.2 uM) did not affect the spiking pattern. 
Unchanged or marginal effects on electrical activity were also seen 
with a specific SK4 channel antagonist, TRAM-34 (1 uM, 
Fig. 1C). However, in some cells TRAM-34 application increased 
the action potential dramatically (Fig. ID) in agreement with 
observations with the SK-channel antagonist apamin [8]. Note 
that before SK-channel block, the cell in Fig. ID was almost 
quiescent, and fired action potentials very infrequently and 
randomly. This increase in spike frequency can be simulated by 
a stochastic version of the model. By including noise in the K(ATP) 
current, an otherwise silent cell produces infrequent action 
potentials evoked by random perturbations (Fig. IE). When the 
SK-conductance is set to 0 nS/pF, the cell starts rapid action 
potential firing driven by the underlying deterministic dynamics. 
The model analysis indicates that this mechanism only works if the 
cell is very near the threshold for electrical activity in the absence 
of the SK-channel antagonist. Dufer et al. [15] suggested a similar, 
important role for SK4 channels in promoting electrical activity in 
murine /J-cells at subthreshold glucose concentrations. Summariz- 
ing, cell-to-cell heterogeneity can explain the differences seen in 
the electrophysiological responses to SK-channel antagonists. 

In addition to spiking electrical activity, human /J-cells often 
show rapid bursting, where clusters of a few action potentials 
(active phases) are separated by hyperpolarized silent phases 
[1,4,9,10,16] (Fig. 2A). The extended model presented here can 
also reproduce this behavior (Fig. 2B) as could the previous version 
of the model [7], where the alternations between silent and active 
phases were controlled by HERG-channels. In contrast, in the 
present version of the model the rapid burst pattern (Fig. 2B, 
upper trace) can be controlled by SK-channels, which in turn are 
regulated by [Ca 2+ ] mem and ultimately by bulk cytosolic Ca 2+ 
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Figure 1. Heterogenous responses to SK-channel block. Note 
the differences in time-scales. A: Simulation, with default parameters, 
showing no effect of SK-channels block (gsK = 0nS/pF during the 
period indicated by the gray bar). B: Experimental recording of spiking 
electrical activity in the same human /J-cell before (left) and during 
(right) application of the SK1-3 channel antagonist UCL-1684 (0.2 \iM). 
C: Experimental recording of spiking electrical activity in the same 



human /J-cell before (left) and during (right) application of the SK4 
channel antagonist TRAM-34 (1 |iM), which had little effect on the 
action potential frequency in this cell. D: Experimental recording of 
spiking electrical activity in the same human /J-cell before (left) and 
during (right) application of the SK4 channel antagonist TRAM-34 
(1 uM), which accelerated the action potential frequency in this cell. E: 
Stochastic simulation reproducing the dramatic effect of SK-channels 
block (gsx = 0 nS/pF during the period indicated by the gray bar). Other 
parameters took default values, except gKATP = 0.0175 nS/pF. 
doi:10.1371/journal.pcbi.1003389.g001 

levels ([Ca 2+ ] c ). The simulated cytosolic Ca 2+ concentration 
shows the characteristic sawtooth pattern (Fig. 2B, lower trace) of a 
slow variable underlying bursting [17,18]. Thus, as in the 
pioneering model by Chay and Keizer [19], [Ca 2+ ] c increases 
during the active phase and activates SK-channels, which 
eventually repolarize the cell. During the silent phase [Ca 2+ ] c 
decreases and SK-channels close, allowing another cycle to occur. 

Na + channels 

Blocking voltage-dependent Na + -channels in human /i-cells 
showing spiking electrical activity with tetrodotoxin (TTX) 
typically reduces the action potential amplitude by —10 mV, 
and broadens its duration [3,9,10] (Fig. 3A). The previous version 
of the model [7] could reproduce these results, though the 
reduction in peak voltage was slightly less than observed 
experimentally. The inclusion of SK-channels in the model leads 
to a greater reduction in the spike amplitude (Fig. 3B, upper trace) 
when Na + -channels are blocked. This improvement is because of 
a mechanism where the slower upstroke in the presence of Na + - 
channel blockers allows submembrane Ca 2+ to build up earlier 
and to higher concentrations (Fig. 3B, lower trace), and 
consequendy to activate more SK-channels, which in turn leads 
to an earlier repolarization reducing the action potential 
amplitude. In other experiments (Fig. 3C) [16], TTX application 
suppresses action potential firing. In agreement, simulated spiking 
electrical activity can be suppressed by TTX application if the cell 
is less excitable because of, for example, smaller Ca 2+ -currents 
(Fig. 3D, upper, black trace). Before TTX application, the 
simulated cell had less hyperpolarized inter-spike membrane 
potential (~ — 61 mV; Fig. 3D) compared to the simulation with 
default parameters (~ — 70 mV, Fig. 3B). This finding is in 
accordance with experimental recordings (compare Fig. 3A and 
3G). The cessation of action potential firing leads to a reduction in 
simulated [Ca 2 + ] mem (Fig. 3D, lower, black trace). The model 
predicts that spiking, electrical activity can continue in presence of 
TTX even in less excitable cells, e.g., with lower depolarizing 
Ca 2+ -currents, if the hyperpolarizing K(ATP)-current is suffi- 
ciendy small (Fig. 3D, upper, gray trace). In this case, [Ca 2+ ] mem is 
nearly unchanged (Fig. 3D, lower, gray trace). Hence, it is the 
relative sizes of the depolarizing and hyperpolarizing currents that 
determine whether TTX application silences the cell or allows the 
cell to remain in a region where action potential firing continues. 
The model thus predicts that in some cells, which stop firing action 
potentials in the presence of TTX, increased glucose concentra- 
tions or sulfonylureas (K(ATP)-channel antagonists) could reintro- 
duce spiking electrical activity. 

More surprisingly, TTX application can change spiking 
electrical activity to rapid bursting in some cells (Fig. 3E). This 
behavior can also be captured by the model (Fig. 3F). To simulate 
this behavior it was necessary to increase the size of the Na + - 
current. Without TTX, the big Na + -current leads to large action 
potentials, which activate sufficient BK-current to send the 
membrane potential back to the hyperpolarized state, allowing a 
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Figure 2. Bursting in human /j-cell. A: Experimental recording of rapid bursting in a human /3-cell. B: Simulation of bursting driven by [Ca 2+ ] 

SK-channels. Default parameters except g s ^ = 0.03 nS/pF, g Kl , = 0.25 nS/pF, nx P Q = -10 mV. 

doi:10.1371/journal.pcbi.1003389.g002 



new action potential to form. With Na + -channels blocked, there 
is insufficient depolarizing current to allow full action potentials to 
develop. In consequence, less BK-current is activated (Fig. 3F, 
lower trace), and the membrane potential enters a regime with 
more complex dynamics where smaller spikes appear in clusters 

from a plateau of 40 mV. The change to bursting activity 

leads to a notable increase in simulated [Ca 2+ ] mem (Fig. 3F, 
middle trace). 



Ca 2+ channels 

High-voltage activated L- and P/Q;type Ca 2+ -currents are 
believed to be directly involved in exocytosis of secretory granules 
in human jS-cells [1,3,4,20,21]. Blocking L-type Ca 2+ -channels 
suppresses electrical activity [3] , which is reproduced by the model 
(Fig. 4A) [7], and the lack of electrical activity is likely the main 
reason for the complete absence of glucose stimulated insulin 
secretion in the presence of L-type Ca 2 + -channel blockers [3] . 
Thus, L-type Ca 2+ -channels participate in the upstroke of action 
potentials and increases excitability of human /J-cells. 

In contrast, and surprisingly, application of the P/GHype Ca 2+ 
-channel antagonist co-agatoxin IVA does not block or slow down 
electrical activity, but leads to an increased spike frequency (Fig. 4B). 
Electrical activity continues also in our model simulations of P/Q- 
type channel block with slightly increased spike frequency (Fig. 4C). 
Reduced Ca 2+ entry leads to lower peak Ca 2+ concentrations in 
the submembrane space ([Ca 2+ ] mem ; Fig. 4D). As a consequence, 
less hyperpolarizing SK-current is activated (Fig. 4E), which leads 
to an increase in spike frequency (Fig. 4G). Hence, the reduction in 
excitability caused by blockage of the P/GMype Ca 2+ -current can 
be overruled by the competing increase in excitability due to the 
smaller SK-current. Experimentally, co-agatoxin IVA application 
reduced the action potential amplitude slightly in 3 of 4 cells (by 



2.0-4.3 mV), a finding that was quantitatively reproduced by the 
model, although the reduction was larger (~7.5 mV in Fig. 4C). A 
direct conclusion from Fig. 4B is that the P/Q_-type Ca 2+ -current is 
not needed for the action potential upstroke, unlike the L-type 
current, probably because of the fact that P/GMype channels 
activate at higher membrane potentials than L-type channels. The 
fact that electrical activity persists with P/Q;type Ca 2+ -channels 
blocked, albeit with lower peak [Ca 2+ ] mem , could underlie the 
finding that (U-agatoxin IVA only partly inhibits insulin secretion [3] . 

Paracrine effects on electrical activity 

The neurotransmitter y-aminobutyric acid (GABA) is secreted 
from pancreatic j8-cells, and has been shown to stimulate electrical 
activity in human /J-cells [12]. In human /i-cells, GABA activates 
GABA^ receptors, which are ligand-gated Cl~ channels, thus 
creating an additional current. Notably, the CI reversal potential 
in human /?-cells is less negative than in many neurons, and positive 
compared to the /?-cell resting potential, which means that CI 
currents, such as the GABA^ receptor current, stimulate action 
potential firing in /?-cells. Hence, GABA is a excitatory transmitter 
in j8-cells, in contrast to its usual inhibitory role in neurons. We 
simulate the addition of GABA by raising the GABA^ receptor 
conductance. In a silent model cell with a rather large K(ATP)- 
conductance, simulated GABA application leads to a single action 

potential whereafter the membrane potential setties at 45 mV 

(Fig. 5A), in close correspondence with the experimental results 
[12]. In an active cell, the simulation of activation of GABA^ 
receptors leads to a minor depolarization and increased action 
potential firing (Fig. 5B), as found experimentally [12]. 

Another neurotransmitter, acetylcholine, might also play a 
paracrine role in human pancreatic islets, where it is released from 
a-cells, and activates muscarinic receptors in /?-cells [22]. 
Muscarinic receptor activation by acetylcholine triggers a 
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Figure 3. Tetrodotoxin (TTX, 0.1 ng/ml) has different effects on electrical activity in human / -cells. A: TTX caused a reduction in action 
potential amplitude in this human /?-cell. B: Simulation with default parameters showing V (upper trace, left axis) and [Ca 2+ ] mem (lower trace, right 
axis), reproducing the data in panel A. C: TTX abolished action potential firing in this human /J-cell. D: Simulations of V and [Ca 2+ ] mcm with default 
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parameters except g Ca L = 0.1 00 nS/pF. With default K(ATP)-channel conductance g K ATP = 0.010 nS/pF, the simulation reproduces the data in panel C 
(black traces). When gKATP = 0.002 nS/pF, the model shows continued firing with Na + -channel block (gray traces). E: TTX changed spiking into rapid 
bursting electrical activity in this human /i-cell. F: Simulation showing V (upper), [Ca 2+ ] mcm (middle), and I BK (lower), reproducing the data in panel E. 
Parameters took default values, except gAr a = 0.7 nS/pF, T„jv„ = 3 ms, g Kv = 0.25 nS/pF, g S K = 0.023 nS/pF, gi eak = 0.012 nS/pF, and nx PQ = —10 mV. 
The extracellular glucose concentration was 6 mM in all experiments. Each couple of experimental traces (panels A, C and E) is from the same human 
/i-cell before (left) and during (right) application of TTX. In the simulations, the Na + -channel conductance g Na was set to 0 nS/pF during the period 
indicated by the gray bars. 
doi:10.1371/journal.pcbi.1003389.g003 



voltage-insensitive Na + -current in mouse pancreatic beta-cells 
[23], and similarly, the muscarinic agonist carbachol activates 
nonselective Na + leak channels (NALCN) in the MIN6 /i-cell line 
[24]. Based on these findings, it was speculated that muscarinic 
activation of NACLN currents in human /J-cells might participate 
in the positive effect of acetylcholine and carbachol on insulin 
secretion [4]. Experimentally, we found that carbachol (20 uM) 
accelerates action potential firing (Fig. 5C). We tested the 
hypothesis of a central role of leak current activation by increasing 
the leak conductance in the model to simulate carbachol 
application, which caused accelerated action potential firing. 
The simulation thus reproduced the experimental data, and lends 
support to the hypothesis that carbachol and acetylcholine can 
accelerate action potential firing via muscarinic receptor-depen- 
dent stimulation of NALCN currents [4]. 

Slow oscillations 

We finally use our model to address the origin of slow rhythmic 
patterns of electrical activity in human j8-cells (Fig. 6A) [4,25], 



which likely underlie slow oscillations in intracellular Ca 2+ [26,27] 
and pulsatile insulin release [28,29]. Based on accumulating 
evidence obtained in rodent islets [5,30], we have previously 
speculated that oscillations in metabolism could drive these 
patterns [7]. In support of this hypothesis, oscillations in ATP 
levels with a period of 3-5 minutes have been observed in human 
/i-cells [13,14]. By adding a glycolytic component [1 1], which can 
oscillate due to positive feedback on the central enzyme 
phosphofructokinase (PFK), our model can indeed simulate such 
periodic modulation of the electrical pattern, where action 
potential firing is interrupted by long silent, hyperpolarized 
periods, which drives slow Ca 2+ oscillations (Fig. 6). 

Discussion 

Human /i-cells show complex and heterogeneous electrophys- 
iological responses to ion channel antagonists. It can therefore 
sometimes be difficult to reach clear conclusions regarding the 
participation of certain ion channels in the various phases of 
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Figure 4. Block of L- and P/Q-type Ca 2+ -channels affects electrical activity differently. A: Spiking electrical activity is suppressed by L-type 
Ca 2+ -channel block in the model with default parameters, and g Ca L = 0 nS/pF during the period indicated by the gray bar. B: Spiking electrical 
activity is accelerated by the application of co-agatoxin IVA in human /i-cells. Recordings from the same human /i-cell in 6 mM extracellular glucose 
before (left) and during (right) application of 200 nM co-agatoxin IVA. C: Model simulation with default parameters of the membrane potential during 
spiking electrical activity under control conditions and after blockage of P/Q-type Ca 2+ -channels (gpg = 0 nS/pF during the period indicated by the 
gray bar). D: In the model, the peak submembrane Ca 2+ -concentration [Ca 2+ ] mem is lower when P/Q-type channels are blocked. E: The reduced 
[Ca 2+ ] mem activate less SK-current when P/Q-type channels are blocked. 
doi:10.1371/journal.pcbi.1003389.g004 
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Figure 5. Paracrine effects on electrical activity. A: Simulation of application of 100 uM GABA to a silent cell (reproducing Fig. 7A in [12]). 
Default parameters except gKATP = 0.021 nS/pF. GABA application was simulated by setting gGABAR to 0.1 nS/pF during the period indicated by the 
gray bar. B: Simulation of application of 10 uM GABA to an active cell (reproducing Fig. 7B in [12]). GABA application was simulated by setting sgabar 
to 0.020 nS/pF during the period indicated by the gray bar. Other parameters took default values. C: Experimental recording of spiking electrical 
activity in the same human /J-cell before (left) and during (right) application of carbachol (20 uM). D: Simulation of accelerated action potential firing 
due to carbachol application. Default parameters except gKATP = 0.016 nS/pF. Carbachol application was simulated by increasing g kak to 0.030 nS/pF 
during the period indicated by the gray bar. 
doi:10.1371/journal.pcbi.1003389.g005 



electrical activity, in particular since some of the electrophysio- 
logical responses are nonintuitive as shown here. A deeper 
understanding of the role of ion channels in electrical activity 
and insulin secretion could have important clinical benefits, since it 
might help in the development of new anti-diabetic drugs. 

We have here shown how mathematical modeling can help in 
interpreting various electrophysiological responses, and in partic- 
ular, to study the effect of competing effects and cell heterogeneity. 
The role of SK-channels in human /?-cells is still not clear. We 
(Fig. 1) and others [8] have found heterogeneous electrophysio- 
logical responses to SK-channel antagonists. Our model suggests 
that these differences can be caused by underlying variations in cell 
excitability: Less excitable /?-cells that produce action potentials 
evoked mostly by stochastic channel dynamics show a clear 
increase in action potential frequency when SK-channels are 
blocked (Fig. IDE). In contrast, spiking electrical activity in very 
active cells is driven by the deterministic dynamics caused by ion 
channel interactions, and is nearly unchanged by SK-channel 
blockers (Fig. 1A-C). We showed also that rapid bursting activity 
can be driven by Ca 2+ and SK-channels (Fig. 2), which could add 
a complementary mechanism to HERG-channel dynamics [7] for 
the control of rapid bursting. 

The wide range of responses to TTX could be accounted for by 
a single model but with different parameters, i.e., differences in the 
relative size of the various currents. A peculiar finding is the 
qualitative change from spiking to rapid bursting seen in some cells 
(Fig. 3E). We suggest that this happens in human /?-cells with large 



Na + -currents. The blockage of this depolarizing current reduces 
the amplitude of the action potentials, and as a consequence, the 
size of the hyperpolarizing BK-current. Under the right condi- 
tions, the combination of these competing events allows the 
membrane potential to enter a bursting regime controlled by SK- 
and/ or HERG-channels (Fig. 3F). Interestingly, it has been found 
that TTX reduces insulin secretion evoked by 6 mM glucose 
greatly, but at glucose levels of 10-20 mM, the effect of TTX on 
secretion is smaller [3,9,10]. Based on our simulations showing 
that less excitable cells cease to fire in the presence of TTX 
(Fig. 3D, black traces), but that lower gKATP can reintroduce 
spiking activity (Fig. 3D, gray traces), we suggest that at low, near- 
threshold glucose levels TTX abolishes electrical activity in many 
cells, which reduces the [Ca 2+ ] mem and consequendy insulin 
secretion greatly (Fig. 3D, black traces). At higher glucose 
concentrations, /?-cells have lower K(ATP)-conductance and in 
some of the cells that stop firing in low glucose concentration the 
effect of TTX on electrical activity and [Ca 2 + ] mem is smaller 
(Fig. 3D, gray traces). Hence, more /i-cells remain active in the 
presence of TTX at high than at low glucose levels. Consequently, 
insulin secretion is more robust to TTX at higher glucose 
concentrations. 

Similarly, insulin release is more affected by the P/ Q;type Ca 2+ 
-channel blocker co-agatoxin IVA at 6 mM (—71%) than at 
20 mM (-31%) glucose [3]. This is in contrast to L-type Ca 2+ - 
channel antagonists, which abolish insulin secretion at both high 
(15-20 mM) and low (6 mM) glucose concentrations [1,3,10]. 
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Figure 6. Metabolically driven slow waves of electrical activity 
and Ca 2+ oscillations. A: Experimental recording of slow oscillations 
in action potential firing in a human /J-cell exposed to 10 mM glucose. 
B-D: Simulation of slow bursting driven by glycolytic oscillations with 
glucose concentration G=10mlvl and default parameters, except 
gift = 0.2 nS/pF, gsK = 0.02 nS/pF, gBK = 0.0l nS/pF. Oscillations in 
glycolysis create pulses of FBP (B), which via ATP production modulates 
K(ATP) channels in a periodic fashion (C). The rhythmic changes in 



K(ATP) conductance drives slow patterns of electrical activity (D), which 
causes oscillations in the intracellular Ca 2+ concentration (E). 
doi:10.1371/journal.pcbi.1003389.g006 

These results concerning L-type Ca 2+ -channel block are easily 
explained by the fact that L-type channel activity is necessary for 
action potential generation [3] (Fig. 4A). In contrast, we showed 
that electrical activity in human /?-cells not only persists, but is 
accelerated by eo-agatoxin IVA (Fig. 4B). The counter-intuitive 
finding of increased excitability and electrical activity when the 
depolarizing P/Q_-type Ca 2+ -current is blocked by eo-agatoxin 
IVA can be accounted for by an even greater reduction in the 
hyperpolarizing SK-current due to reduced Ca 2 + -influx and 
consequently lower [Ca 2+ ] mcm . 

Our mathematical modeling confirmed that GABA released 
from human /?-cells can have a role as a positive feedback 
messenger. GABA application has been shown to depolarize both 
silent and active human /?-cells [12], which was reproduced here. 
A detailed characterization of GABA^ receptor currents would 
refine the analysis presented here. 

Data from mouse /?-cells [23] and the MIN-6 jS-cell line [24] 
suggest that muscaric agonists such as carbachol and acetylcholine 
stimulate insulin secretion partly by activating NACLN currents. 
Using our model we could translate this finding to the human 
scenario, thus testing the hypothesis that this mechanism is also 
operating in human /J-cells [4]. Our simulations confirmed that 
increased leak currents can underlie the change in electrical 
activity found experimentally (Fig. 5C). The incretin hormone 
glucagon-like peptide 1 (GLP-1) has also been shown to act pardy 
via activation of leak channels [31], a mechanism which might be 
involved in activating otherwise silent /?-cells [7,32,33]. These 
results suggest that leak currents could play important roles in 
controlling electrical activity in /?-cells, and potentially be 
pharmacological targets. Further studies are clearly needed to 
investigate these questions. 

We were also able to simulate slow rhythmic electrical activity 
patterns by adding an oscillatory glycolytic component to the 
model. To date, there is to our knowledge no evidence of 
oscillations in glycolytic variables in human (or rodent) /?-cells or 
islets, but ATP levels have been found to fluctuate rhythmically 
also in human /?-cells [13,14], supporting the idea of metabolism 
having a pacemaker role. In agreement, data from rodent /8-cells 
show accumulating evidence for oscillations in metabolism playing 
an important role in controlling pulsatile insulin secretion [5,30]. It 
will be interesting to see if these findings in rodents are applicable 
to human /J-cells. 

Regarding the model development, the inclusion of SK- 
channels in the model provided insight that was not within reach 
with the previous version of the model [7]. Besides the direct 
investigation of the role of SK-channels, the acceleration in action 
potential firing seen with P/ Ojtype Ca 2 + channel blockers (Fig. 4) 
can not be reproduced by the older version of the model without 
SK-channels [7]. Moreover, considering the effect of TTX on 
spike amplitude, a better correspondence between experiments 
and simulations was found with SK-channels included in the 
model. To model SK-channel activation accurately, we made a 
special effort to describe [Ca 2+ ] mem carefully. Submembrane 
Ca 2+ responds rapidly to an action potential, while [Ca 2 + ] c 
integrates many action potentials. The rapid submembrane 
dynamics has important consequences for the study of the role 
of SK-channels in spiking electrical activity, e.g., it was crucial for 
explaining the larger effect of TTX on spike amplitude in this 
version of the model. Most models of electrical activity in rodent /?- 
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cells do not include a submembrane Ca 2+ compartment, but these 
models were typically built to explain the slow bursting patterns 
seen in rodent islets with a period of tens of seconds. For these long 
time scales, the rapid dynamics in the submembrane compartment 
is not important. In contrast, the situation is different in human /?- 
cells with their faster dynamics. 

Methods 

Modeling 

We build on the previously published Hodgkin-Huxley type 
model for human /?-cells [7], which was mainly based on the 
results of Braun et al. [3], who carefully assured that investigated 
human islet cells were jS-cells. We include SK-channels in the 
model. Since these channels are Ca 2+ -sensitive and located at 
some distance from Ca 2 + -channels [34] we also model Ca 2+ - 
dynamics in a submembrane layer controlling SK-channel 
activity. 

The membrane potential V (measured in mV) develops in time 
(measured in ms) according to 



dV 



— — iJsK + Ibk + Ikv + Iherg + Inu 



dt (1) 

+ IcaL + IcaPQ + IcaT + IrA TP + 1 leak + IgABAr) ■ 

All currents (measured in pA/pF), except the SK-current Isk and 
the GABA^ receptor mediated current Igabar, are modeled as in 
[7]. Expressions and parameters are given below. For the 
stochastic simulation in Fig. IE, we included "conductance noise" 
[35] in the K(ATP) current by multiplying Ikatp by a stochastic 
factor (l+0.2r ( ), where T t is a standard Gaussian white-noise 
process with zero mean and mean square <r ( ,r. s > = S(t — s), see 
also [36-38]. 

SK-channels are assumed to activate instantaneously in 
response to Ca 2+ elevations at the plasma membrane but away 
from Ca 2+ channels [34], and are modeled as [39] 



IsK=gSK 



CdL 



-XV-Vk). 



(2) 



In human jS-cells, flash-released Ca 2+ triggered a ~10 pA current at 
a holding current of — 60 mV, presumably through SK-channels [8] . 
Assuming that SK-channels were nearly saturated by Ca 2+ , the 
maximal SK-conductance is estimated to be gs/f* 10pA/ 
(-60mV- V K )/C m wO.l nS/pF. Here, C m = 10 pF is the capac- 
itance of the plasma membrane [3] . 

In Eq. 2, Ca m is the submembrane Ca 2+ concentration 
([Ca 2+ ] mem ; measured in uM), which is described by a single 
compartment model [21] 



=faC m ( — IcaL — IcaPQ —IcaT)/ Vol m 



dCa„ 

~dT " ■"- — ™ • (3) 

-f{Vol c /Vol m )[B{Ca m - Ca c ) + (JpMCA +Jncx% 

where / = 0.01 is the ratio of free-to-total Ca 2+ , 
a = 5.18 x 10~ 15 umol/pA/ms changes current to flux, and 
Vol m and Vol c are the volumes of the submembrane compartment 
and the bulk cytosol, respectively. B describes the flux of Ca 2+ 
from the submembrane compartment to the bulk cytosol, Jrmca is 
the flux through plasma membrane Ca 2+ -ATPases, and Jncx 
represents Ca 2+ flux through the Na + - Ca 2+ exchanger. 
Cytosolic Ca 2+ (Ca c ; measured in uM) follows 



dCa c 
dt 



=f[B{Ca m — Ca c ) — JsERCA+Jleak], 



(4) 



where Jserca describes SERCA pump-dependent sequestration 
of Ca 2+ into the endoplasmic reticulum (ER), and Jy a k is a leak 
flux from the ER to the cytosol. Expressions and parameters for 
the Ca 2+ fluxes are taken from [40]. 

The submembrane compartment volume is estimated based on 
the considerations of Klingauf and Neher [41], who found that a 
shell model (in contrast to a domain model) describes submem- 
brane Ca 2+ satisfactorily when the shell-depth is chosen correctly. 
The Ca 2+ dynamics between channels can be estimated from a 
shell model at a depth of ~23% of the distance to a Ca 2+ - 
channel. In mouse /i-cells the interchannel distance has been 
estimated to be ~1200 nm [42]. Moreover, SK-channels are 
located >50 nm from Ca 2+ channels [34]. 

Based on these considerations, we modeled the submembrane 
space controlling SK-channels as a shell of depth ~ 190 nm. The 
radius of a human j3-cell is ~13 um, which gives cell volume 
{Vole), shell volume (Vol m ) and internal surface area (A m ) of the 
shell, of 

Ko/„ = 1.15pL=1150um 3 , Ko/ m = 0.1pL, ^l m = 530um 2 .(5) 
The flux-constant B can then be calculated as [43] 

A m 



B = D ( 



VoLd m 



(6) 



where d m is a typical length scale. We set d m to 1 u.m, which 
together with the diffusion constant for Ca 2+ , Z) Cfl = 220u.m 2 /s 
[41,44], gives 5 = 0.1 ms- 1 . 

In human j8-cells, GABA activates GABA^ receptors, which are 
ligand-gated CI channels. We model the current carried by 
GABA/( receptor as a passive current with the expression 



Igabar— gGABAp{V — Vci), 



(7) 



where gGABAR is the GABA^ receptor conductance, and 
Vci = — 40 mV is the chloride reversal potential [4] . We estimate 
gGABAR from the findings that 1 mM GABA evokes a current of 
9.4 pA/pF (but with substantial cell-to-cell variation) at a holding 
potential of —70 mV [12], which yields a conductance of 
~0.3 nS/pF. To simulate the changes in firing patterns evoked 
by lower GABA concentrations (10 or 100 uM) [12], we take into 
consideration the does-response curve [45] for the al p3 y2 
subunits, which are the most highly expressed subunits in human 
p-celh [12]. At 10 uM the GABA-evoked current is > 10-fold 
smaller compared to 1 mM GABA, and we set gGABAR = 0.02 nS/ 
pF. At 100 uM, the reduction is about 2-fold compared to 1 mM. 
We used gGABAR =0.10 nS/pF to simulate application of 100 uM 
GABA. 

To investigate slow electrical patterns (Fig. 6) we added a 
glycolytic component [1 1], which drives ATP levels and K(ATP) 
channel activity. The glycolytic subsystem can oscillate due to 
positive feedback on the enzyme phosphofructokinase (PFK) from 
its product fructose- 1,6-bisphosphate (FBP). The glycolytic equa- 
tions are 



d G6P F6P 
dt 



VgK — VpFK, 



(8) 
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clFBP 
dt 



= V P 



Vfba , 



TRAM- 3 4 from Sigma-Aldrich (Oakville, ON, Canada). Figures 
with experimental responses to ion channel antagonists (Figs. 1, 3, 
(9) 4 and 5) show recordings from the same cell before (ctrl) and after 
application of the blocker. 



d DHAPG3P 
dt 



■2V FB a-V g 



(10) 



where Vgk is the rate of glucokinase, which phosphorylates 
glucose to glucose-6-phosphate (G6P). G6P is assumed to be in 
equilibrium with fructose-6-phosphate (F6P), the substrate for 
PFK, and G6PF6P is the sum of G6P and F6P. V PFK is the rate 
of PFK producing FBP, which is subsequently removed by 
fructose-bisphosphate aldolase (FBA), which produces glyceralde- 
hyde-3-phosphate (G3P) and dihydroxyacetone-phosphate 
(DHAP) with rate Vfba- DHAP and G3P are assumed to be in 
equilibrium, and DHAPG3P indicates their sum. Finally, G3P 
serves as substrate for glyceraldehyde-3-phosphate dehydrogenase 
(GAPDH with rate Vqapdh), which via the lower part of glycolysis 
eventually stimulates mitochondrial ATP production. We intro- 
duce a phenomenological variable a that mimics ATP levels, and 
is model by 



Model equations and parameters 

For completeness, we report all expressions and parameters of 
the mathematical model here. For details, please refer to the 
Modeling section above and the previous article [7]. 

The main variables, membrane potential, V, submembrane 
Ca 2+ , Ca m , and cytosolic Ca 2+ , Ca c , are described by 



dV 



= — (Isk + Ibk + Ikv + Iherg + iNa 



dt (13) 

+ IcaL + IcaPQ + IcaT + IkATP + heak + IqABAr)- 



dCa„ 



=faC m ( - I C aL - IcaPQ ~ IcuT)/ Vol m 



(14) 



-/( Void Vol m )[B(Ca„, - Ca c ) + (J PM ca + Jncx)] , 



da 

777 — * GAPDH— K-A «■ 

dt 



(11) 



The K(ATP) conductance depends inversely on a, and is modeled 

as 



Skatp = &KATP / (\ + a). 



(12) 



Expressions and parameters are given below. 

Simulations were done in XPPAUT [46] with the cvode solver, 
except the stochastic simulation in Fig. IE, which was performed 
with the implicit backward Euler method. Computer code can be 
found as supplementary material, or downloaded from http:// 
www.dei.unipd.it/ pedersen. 

Experiments 

Human pancreatic islets were obtained with ethical approval 
and clinical consent from non-diabetic organ donors. All studies 
were approved by the Human Research Ethics Board at the 
University of Alberta. The islets were dispersed into single cells by 
incubation in Ca 2+ free buffer and plated onto 35 mm plastic 
Petri dishes. The cells were incubated in RPMI 1640 culture 
medium containing 7.5 mM glucose for >24 h prior to the 
experiments. Patch-pipettes were pulled from borosilicate glass to 
a tip resistance of 6-9 Mfi when filled with intracellular solution. 
The membrane potential was measured in the perforated-patch 
whole-cell configuration, using an EPC-10 amplifier and Patch- 
master software (HEFvA, Lambrecht, Germany). The cells were 
constantly perifused with heated bath solution during the 
experiment to maintain a temperature of 31— 33°C. The 
extracellular solution consisted of (in mM) 140 NaCl, 3.6 KC1, 
0.5 MgSOz,, 1.5 CaCl 2 , 10 HEPES, 0.5 NaH 2 P0 4 , 5 NaHC0 3 
and 6 glucose (pH was adjusted to 7.4 with NaOH). The pipette 
solution contained (in mM) 76 K 2 S0 4 , 10 KC1, 10 NaCl, 1 
MgCl 2 , 5 HEPES (pH 7.35 with KOH) and 0.24 mg/ml 
amphotericin B. /?-cells were identified by immunostaining (18 
out of 28 cells) or by size when immunostaining was not possible 
(cell capacitance >6 pF, [3]). Tetrodotoxin (TTX) and co-agatoxin 
IVA were purchased from Alomone Labs (Jerusalem, Israel), 
UCL-1684 was obtained from R&D Systems (Minneapolis, MN), 



dCa c 
dt 



=f[B(Ca m - Ca c ) - Jserca + J leak] (15) 



The currents are 



Ca" 



SK ~ ^ u m 



(V-V K ), 



IcaPQ = gCaPQmCaPQ,oo( V)( V - V Ca ), 



IcaT=gCaTm Ca T,aj(V)h C aT(V- V C a), 



IlC(ATP) = gK(ATP){V — Vk), 



heak = gteaki V — Vl ea k), 



IGABAR — gGABARi V — Va), 



(16) 



Ibk = gBKm BK ( - I Ca ( V) + B BK ){ V— Vk), (17) 



lKv=gKvm Kv (V-V K ), (18) 



Iherg = gHERG^HERG^HERGi V — Vk), (19) 



iNa = gNam Na ,a,( V)ll Na { V - V Na ), (20) 



ICaL = gCaLm C aL,oo( V)h C aL.{ V ~ V C a), (21) 



(22) 



(23) 



(24) 



(25) 



(26) 
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Table 1. Default model parameters. 



Parameter 






Ref. 


Parameter 






Ref. 


V K 


-75 


mV 


[3] 


v Na 


70 


mV 


[7] 


Vca 


65 


mV 


[7] 


V a 


-40 


mV 


[4] 


9sk 


0.1 


nS/pF 


[8] 




0.57 


iiM 


[39] 


n 


5.2 




[39] 










&BK 


0.020 


nS/pA 


[3] 


TmBK 


2 


ms 


[3] 




0 


mV 


[3] 


n mBK 


-10 


mV 


[3] 


&BK 


20 


pA/pF 


[3] 










9kv 


1.000 


nS/pF 


[3] 


T mKv,0 


2 


ms 


[3] 


VmKv 


0 


mV 


[3] 


n mKv 


-10 


mV 


[3] 


9herg 


0 


nS/pF 


+ [7,47] 










V/nHERB 


-30 


mV 


[47] 


n mHERG 


-10 


mV 


[47] 


VhHERG 


-42 


mV 


[47] 


n hHERG 


17.5 


mV 


[47] 


z mHERG 


100 


ms 


[48] 


T hHERG 


50 


ms 


[47] 


9Na 


0.400 


nS/pF 


[3] 


T hNa 


2 


ms 


[3] 


VmNa 


-18 


mV 


[3] 


n mNa 


-5 


mV 


[3] 


VhNA 


-42 


mV 


[3] 


n hNa 


6 


mV 


[3] 


9caL 


0.140 


nS/pF 


[3] 


thCaL 


20 


ms 


[7] 


VmCaL 


-25 


mV 


[3] 


n mCaL 


-6 


mV 


[3] 


9caPQ 


0.170 


nS/pF 


[3] 










VmCaPQ 


-10 


mV 


[3] 


n mCaPQ 


-6 


mV 


[3] 


9caT 


0.050 


nS/pF 


[3] 


thCaT 


7 


ms 


[3] 


VmCaT 


-40 


mV 


[3] 


^mCaT 


-4 


mV 


[3] 


VhCaT 


-64 


mV 


[3] 


n hCaT 


8 


mV 


[3] 


9k(atp) 


0.010 


nS/pF 


[1] 


Qgabar 


0 


nS/pF 


[12,45] 


9leak 


0.015 


nS/pF 


[1] 


Vleak 


-30 


mV 


[7] 


J SERCA.max 


0.060 


^iM/ms 


+ [40] 


K-SERCA 


0.27 


(iM 


[40] 


JpMCA.max 


0.021 


(iM/ms 


[40] 


KpMCA 


0.50 


HM 


[40] 




0.00094 


jiM/ms 


[40,49] 


Jncx.o 


0.01867 


ms~' 


[40,49] 


f 


0.01 






Vol c 


1.15 


x10~ 12 L 




B 


0.1 


ms" 1 




Vol m 


0.1 


x10~ 12 L 




a 


5.18x10~ 15 jimol/pA/ms 
















0.0000556 


mM/ms 


[11] 




8 


mM 


[11] 


r>GK 


1.7 




[11] 


G 


10 


mM 


+ 




0.000556 


mM/ms 


[11] 


KpFK 


4.0 


mM 


[11] 


"PFK 


2.5 






hact 


1 




[11] 


XpFK 


0.01 


mM 




h x 


2.5 




[11] 


a-G 


5.0 














VfBA.mhx 


0.000139 


mM/ms 






0.005 


mM 


[11] 


P FBA 


0.5 


mM 




Ofba 


0.275 


mM 


[11] 


^GADPH.max 


0.00139 


mM/ms 




Kgadph 


0.005 


mM 


[11] 
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Table 1. Cont. 







Parameter 






Ref. 


Parameter 






Ref. 




0.3 




[11] 


Km 


0.045455 




[11] 




0.0001 


ms-' 


+ 


gKATP 


0.050 


nS/pF 


+ 



Default model parameters used in the manuscript unless mentioned otherwise. Parameter values are based on the indicated literature references (+ indicates adjusted 
parameters), see also reference [7] for a discussion of the parameters introduced there. All glycolytic parameters are taken without modification from reference [11], 
where a discussion of their values based on experimental data can be found. HERG channel conductance, gHERd was set to zero in the present work to investigate 
whether SK-channels can substitute for HERG channels, e.g., in driving bursting. Based on [47], the previous version of the model [7] used ghERG = 0.2 nS/pF. The 
conclusions presented here are not sensitive to whether or not HERG currents are included. 
doi:10.1371/journal.pcbi.1003389.t001 



where 



ICa( V) = IcaL + ICaPQ + ICaT, 



(27) 



d G6PF6P 



VqK — VpFK, 



(35) 



and activation variables (and similarly inactivations variabels, hx, 
where X denotes the type of current) follow 



dm x = m x , m {V)-m x 
dt x mX 



(28) 



where x m x (respectively ihx) is the time-constant of activation 
(respectively inactivation for hx), and Wx,oo(f0 (respectively 
hx,co(V)) is the steady-state voltage-dependent activation (respec- 
tively inactivation) of the current. The steady-state activation (and 
inactivation) functions are described with Boltzmann functions, 



m x ^(V)-- 



1 

l + exp((V-V m x)/n m xY 



except 



hcaL,aj( V ) = 

max(0,min{ 1,1 + [m Cu LM V){ V - Fo»)]/57 mV}), 



(29) 



(30) 



for Ca 2+ -dependent inactivation of L-type Ca 2+ channels. The 
time-constant for activation of Kv-channes is assumed to be 
voltage-dependent [3,7], 



' W v ,o+10exp( 2 ^ m y, )ms, for V>26.6mV 
x mK y=\ \ 6mV I (31) 



\tmKv,o + 30 ms, 



Ca 2+ -fluxes are [40] 



for V< 26.6 mV. 



d FBP 



d DHAP-GiP 
dt 



Vfba, 



PFK y FBA 



'■ 2 Vfba — Vgapdh, 



(36) 



(37) 



which controls the electrophysiological subsystem via the "ATP- 
mimetic" a and K(ATP)-channels, as described by 



Here 



da ; 

-77 — V GAPDH — HA « , 

dt 



gKATP = gKATP I '(1 +fl). 



V G k=V, 



GK— * GK.max 



K'g, K + Gi< GK ' 



(38) 



(39) 



(40) 



Vpfk = VpFKmia 



(F6P\" (FBP} 

\KpFKj 



( FBP\ 

(F6Py FBP > } 1+ {xpfk) 
\Kpfk) (FBP\ hx h( F BP) 



(41) 



J SERCA — J SERCA,mzx ~js2T 



Cat 



■CaV 



CCIn 



'PMCA,m-dx 



K P 



i + Ca„ 



(32) 



(33) 



Vfba. 



WAX I 



FBP G3P x DHAP 



\Kfba PfbaQfbaK-fba 

FBP _|_ DHAP _j_ G3PxDHAP 
K FBA QFBA P FBA QpBA 



(42) 



Jncx — Jncx.o Ca„ 



(34) 



VgADPH = VaADPH.m-M 



G3P 



K G 



G3P' 



(43) 



Glycolysis was modeled by [1 1] 



where 
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F6P=(G6PF6P)K GPI /(l+K GPI ), (44) 



G3P=(DHAPG3P)K T pi/(l+K T pi), (45) 



DHAP = {DHAPG3P)-G3P, (46) 

and 



h(FBP) = h PFK -(h PFK -h act )- (47) 

J^FBA -\-rJsr 

Default parameters are given in Table 1. 
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